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Abstract 



We introduce a three-dimensional (3D) model of optical media with the quadratic (x^^'') nonlin- 
earity and an effective 2D isotropic harmonic-oscillator (HO) potential. While it is well known that 
3D x'^' solitons with embedded vorticity ("vortical light bullets") are unstable in the free space, we 
demonstrate that they have a broad stability region in the present model, being supported by the 
HO potential against the splitting instability. The shape of the vortical solitons may be accurately 
predicted by the variational approximation (VA). They exist above a threshold value of the total en- 
ergy (norm) and below another critical value, which determines a stability boundary. The existence 
' threshold vanishes is a part of the parameter space, depending on the mismatch parameter, which is 

explained by means of the comparison with the 2D counterpart of the system. Above the stability 
boundary, the vortex features shape oscillations, periodically breaking its axisymmetric form and 
restoring it. Collisions between vortices moving in the longitudinal direction are studied too. The 
collision is strongly inelastic at relatively small values of the velocities, breaking the two colliding 
C/3 ■ vortices into three, with the same vorticity. The results suggest a possibility of the creation of stable 

^ ' 3D optical solitons with the intrinsic vorticity. 

Oh' 



I. INTRODUCTION 



The second-harmonic-generating {x^^\ quadratic) nonlinearity is one of basic features of photonic media. Among 
other fundamental effects, the x^^^ interactions give rise to two-color solitons of diverse types In particular, 

the use of this nonlinearity is critically important for the creation of stable two- and three-dimensional (2D and 3D) 
solitons, because, in contrast to the cubic (x^'^-*) self-focusing, the x^^-* interaction between the fundamental-frequency 
(FF) and second-harmonic (SH) fields does not induce the wave collapse [^, which destabilizes multidimensional 
Tij" , solitons in x''^'* settings 0, (7|. Due to this circumstance, the x*'^'' solitons were originally created, in the spatial 
• domain, as stable (2-|-l)-dimensional beams Q. The possibility of the creation of fully localized 3D spatiotemporal 
' solitons is suggested too by the absence of the collapse @-0- This has not been achieved yet in experiments, the best 
result being a spatiotemporal soliton which is self-trapped in the longitudinal and one transverse directions, while 
T-H ^ the confinement along the other transverse coordinate was provided by the waveguiding structure (isl. [l^. The x'^'* 
i nonlinearity is also promising for the creation of ultrashort few-cycle optical solitons [15| . 

The creation of self-trapped beams with intrinsic vorticity is another natural possibility in the 2D setting, with 
the SH and FF fields carrying topological charges m and m/2, respectively. For even m, these modes are classified 
as solitary vortices with topological charge m/2. States with odd values of m are not possible, as the topological 
charge of the FF component cannot take half-integer values (vortices with a half-integer optical angular momentum, 
in the form of mixed screw-edge dislocations, can be created by passing the holding beam through a spiral-phase 
plate displaced off the beam's axis Unlike the fundamental x*''^'' solitons with m = 0, solitary vortices in the 

free space are always unstabl e ag ainst azimuthal perturbations, which split them into sets of separating segments, as 
shown both theoretically p7H2Cll | and experimentally [2l|. A similar splitting instability of vortices was predicted in 
the framework of the three- wave (alias Type-II) x^^'' system, which includes two components of the FF field [23|. 

Vortical solitons are well known too as solutions to the 2D and 3D nonlinear Schrodinger (NLS) equation with the 
cubic (x^'^^) self- focusing nonlinearity They are also unstable against azimuthal perturbations, which develops 
faster than the above-mentioned collapse, i.e., the vortex splits into fragments, which later suffer the intrinsic collapse 
Q . A method for the stabilization of the vortices was elaborated in the framework of the 2D and 3D Gross-Pitaevskii 
(GP) equations for atomic Bose-Einstein condensates (BECs) with attractive inter-atomic interactions, which are 
similar to the NLS equations for photonic media with the self- focusing x^'^'' nonlinearity [1^. It has been shown 
in detail theoretically [26n32i] that the axisymmetric harmonic-oscillator (HO) potential stabilizes the fundamental 



solitons in their whole existence domain, and vortices with topological charge 1 in a part of the region where they 
exist. 

In optical media, the effective 2D trapping potential can be realized too, in the form of an appropriate transverse 
profile of the local refractive index [s*] . The stabilization of 2D vortex solitons by means of the isotropic HO potential 
in the medium with the x^^^ nonlincarity was reported in our recent work .33.|. A 2D photonic crystal can be used to 
build this setting, with the x^^^ nonlincarity provided by the poled liquid jSJ] or solid [35[ material filling its voids. In 
fact, the profile of the effective potential does not need to be exactly parabolic (HO), due to the strong localization 
of the trapped modes. Our recent analysis _33] demonstrates the stabilization of both two-color x^"^^ vortices and 
single-color ones, only the SH component being present in the latter case. Stability regions were found for both 
types of the modes. In particular, while the single-color states are actually linear modes supported by the trapping 
potential, their stability is determined by the x^^-* interactions with small perturbations in the FF field. Two-color 
trapped vortices exist above a finite threshold values, which is identical to the instability boundary of the single-color 
vortex mode. 

A natural extension of the previous analysis, which is the subject of the present work, is to consider a possibility of 
the stabilization of two-color 3D vortex solitons ( "vortical light bullets" 3] ) by the 2D axisymmetric trapping potential, 
which may open the way to creating such self-trapped spatiotemporal modes. Thus far, they evaded attempts of the 
experimental observation (it is obvious that, unlike the 2D setting, single-color localized modes cannot exist in the 
3D x^^^ system with the 2D potential). The 3D model, based on the system of coupled equations for the FF and 
SH fields, is introduced in Section II. In particular, an essential parameter, specific to the 3D model, is the ratio of 
the group-velocity-dispersion (GVD) coefficients for the FF and SH waves. In fact, essentially the same system of 
GP equations for atomic and molecular wave functions, applies to the BEG in atomic-molecular mixtures, where the 
role of the GVD ratio is played by the ratio of the atomic and molecular masses [3^ - [39j , hence the mechanism of the 
stabilization of the 3D vortex solitons trapped in the 2D HO potential may be implemented in the BEG mixture too. 
In Section II we also formulate the variational approximation (VA) for stationary states. 

Stationary solutions for the trapped 3D vortical modes are considered in Section HI. The results are obtained by 
means of numerical methods, and in a semi-analytical form based on the VA. In those cases when the underlying 
ansatz is appropriate, the accuracy of the VA is very good, in comparison with numerical results. Depending on the 
mismatch parameter, the vortex modes may exist above a finite threshold value of the norm (total energy), which can 
be explained in a simple form. The stability of the vortex solitons is considered, by means of systematic simulations 
of the perturbed evolution, in Section IV. The vortex modes are stable below a certain critical value of the norm. 
Above the instability boundary, the vortex develop oscillations, maintaining its vorticity and demonstrating periodic 
breaking and recovery of the axial symmetry of its shape. Gollisions between vortex solitons moving in the longitudinal 
direction are also considered in Section IV, a noteworthy result being that, at sufficiently small velocities, two colliding 
vortices may break up into three. The paper is concluded by Section V. 



II. THE MODEL AND THE VARIATIONAL APPROXIMATION 



The equations for the evolution of the FF and SH field amplitudes, u (x, y, z, r) and v (x, ?/, 2, r), in the 3D setting, 
with propagation distance z, transverse coordinates (x, ?/), and reduced time t = t — z/V are taken in the known form 
[ll-Q {y is the common group velocity of the FF and SH waves, assuming that the group-velocity mismatch between 
them may be made equal to zero by means of an appropriate Galilean boost): 

.du 1 / 32 32 92 X , , 

2^^+2(;^ + ^ + ^^j-'^- + r'-4t/(x,.). = 0, (I) 

where the asterisk stands for the complex conjugate, q is the real mismatch coefficient, the anomalous-GVD coefficient 
at the FF is normalized to be 1, and D is the ratio of the SH and FF GVD coefficients. In the case of Z? = I, the 
diffraction-dispersion operators for the FF and SH components feature the spatiotemporal isotropy 0. As said 
above, the axisymmetric modulation of the refractive index is represented by a spatially-isotropic HO potential, 
U(x,y) = ($7^/2) (a;^ +y^)- By means of a further rescaling, we fix = 1/2, while q is kept as a free constant. 
Stationary modes are characterized by their total energy (norm), 

I I I {\uf + 4:\vf) dxdydr, (2) 



which is dynamical invariant of the system. Three other invariants are the angular momentum: 

M = J J J Im{u* {yux — xuy)} dxdydT, 
the linear momentum in the t direction: 



P = yyy lm{u*ur}dxdyd7 



(3) 



(4) 



and the Hamiltonian: 



H = 



(I 



Map + lUyP + |MtP + \Vx 



\Vy\^ +D\Vr 



+U{r) + 4|wp) + q\v\'^ ~ ^ {u'^v* + u^*v) | dxdydr. 
Equations ([T]) can be derived from the corresponding action, S = J Ldz, with Lagrangian 



(5) 



L 



-U{r) + 4|f |2) - q\v\^ + - {u^v* + u^*v) \ dxdydr. 



(6) 



This fact may be used to develop the VA for spatiotemporal vortical solitons with propagation constant fc, based on 
the following natural ansatz, which corresponds to topological charge 1 and is written in terms of polar coordinates 
(r, 6) in the plane of {x,y): 



uqt exp (ikz — air^ + iO) vqt^ exp (ikz — 02^^ + 2i0) 
u ~ — . , „ , -, V = 



cosh (/3t) 



cosh {(3t) 



(7) 



Here (uo,vo) and ai are, respectively, the amplitudes and radial widths of the two components, and P ^ is their 
common temporal width. The norm ^ of the ansatz is 

N = 27r dr H rdr{\uf + A\vf) ^^(^ + ^) (8) 



The substitution of ansatz ^ into the Lagrangian yields 



L 
2^ 



kN 



ulP DvlP 



2ti 2ai/3 4a2/3 24a? 24ai 
"8ap ~ Aa\P ~ 4^ ^ 2l3{2ai + a2f' 



(9) 



Using the relation following from Eq. ([8]), Vq ~ [A^o/? — Uo/(4<^i)] "Ij '^ith A^o = N/{2tt), Lagrangian dH) is cast 
into the form in which A^o raay be considered as a given constant, while uo,ai,a2, and /3 are treated as variational 
parameters: 



2n 2aiP 4/3 V ^(4) 24ai 24 



(■ 



2 3/2 

^0"2 



4a? y 8a?/3 
2 \ 1/2 



(10) 



V4a2/3 ' 4/?; V"'''' 4a? y ' 2/3(2ai + a2)3 
The system of variational equations produced by the latter Lagrangian, dL/duo = dL/dai = dL/da2 ~ dL/d/3 = 



take a cumbersome form, which, nevertheless, admits a numerical solution: 

3uoa2 uo/3 D/Suq rt^uo Sfi^uo 




III. NUMERICAL AND SEMI- ANALYTICAL RESULTS FOR STATIONARY VORTICES 

At first, we numerically solved Eqs. ([1]) by means of the imaginary-time-evolution method (it was earlier applied to 
searching for vortical trapped states in the framework of the GP equation with the cubic nonlinearity (40|'). starting 
with initial conditions corresponding to the 3D vortical wave form with topological charge 1 such as 

uix, 2/, r) = 0.5re~"-25'^%-"-5^'e*^ v{x, y, r) = 0. (12) 

Shapes of the so generated localized vortical states were compared to the prediction of the VA, obtained for the same 
parameters, and the stability of the states was checked by simulations of their perturbed evolution in the framework 
of Eq. H]) in real time. Equations ([T]) were solved numerically in the 3D domain defined as < {x, y, r} < {L, L, T} , 
with LxLxT = 25 X 25 X 50, and periodic boundary conditions. Accordingly, the polar coordinates (r, 6) are defined 
with respect to the central point, x = y = L/2. This size of the integration domain was sufficient to completely 
accommodate all the vortical solitons produced by the numerical solution. 

A characteristic example is presented in Fig. [1] which shows a plot of \u{x,y,T)\ in cross sections drawn through 
T = Lr/2 (a) and y = L/2 (b) at = 40, q = 1, D = l,n = 0.5. In particular. Fig. [TJa) shows the vortex structure per 
se, with amplitude |u| vanishing at the center, and Fig. [Ijb) shows the self-trapped structure in the r direction. For the 
same parameters, the numerical solution of the variational equations (ITT|) yields uq = 0.704, ai = 0.268, a2 — 0.507, 
and /3 = 0.438. To compare the full numerical results with the prediction of the VA, the solid curve in Fig.[2][a) shows 
the numerically found profile of \u{x ~ L/2, y,T)\ aX y = L/2, r = Lr/2, and the dashed curve shows the VA-predicted 
counterpart of this profile, taken as per Eq. JT]), i.e., uq\x\ exp (— aix^). Further, the solid curve in Fig. ^b) shows 
the temporal profile, \u{x, y,T~ Lr/2)\, at fixed x — 1.663 and y — L/2, and the dashed curve shows its VA-predicted 
counterpart, uo|a;| exp (— aix^) /cosh(/3r) at the same point, x — 1.663. It is seen that the VA predicts the numerical 
results in a practically exact form. 

Systematic results obtained for the vortex-soliton family in the direct numerical and VA forms are summarized in 
Fig. [3l in which panel (a) shows the slope of the vortical profile at its center, 

slope = r~'^\u{r,T)\ at r — ;> 0, (13) 

as a function of the total energy (norm) N , for two values of the mismatch, q = 1 and 0.6. In the latter case, the 
vortex soliton disappears discontinuously at a critical value N = Nd ~ 19.5, as demonstrated by the numerical results 
and the VA alike, and the vortices do not exist at TV < A^ci- Further, Fig. ^h) shows the critical value as a function 
of the mismatch. Note that Nd vanishes at q > 1. 




FIG. 1: Plots of \u{x,y,T)\ in the cross section drawn through r — T/2 = 25 (a) and y = 12.5 (b). 



(a) 



(b) 





FIG. 2: (a) The plot of \u{x — L/2,j/,r)| at y = L/2,t = T/2 obtained from the direct numerical solution (the solid curve) 
and its VA-predicted counterpart, iio|a;| exp (— Qia; ), with Mo = 0.704, ai = 0.268 (the dashed curve), (b) The same for 
\u{x,y,T — T/2)\ at X — 1.663, y — L/2, the respective VA-predicted profile being uolxl exp (— aix^) sech(/3r) with ug = 
0.704, Qi = 0.268, /3 = 0.438. 



It is relevant to note that a finite threshold value of the norm, N — N^^ ' , necessary for the existence of two-color 
vortical solitons in the 2D version of Eqs. ([T]), was found in our recent work [s^. As well as in the present model, 
^ci^'' vanishes at 5 > 1. The difference from the 3D case is that stationary 2D single-color vortices trapped in the 
HO potential, with the zero FF component (i.e., these are linear SH modes), exist at iV < Nd- Obviously, this is 
not possible in the 3D case, as the 2D potential cannot trap 3D linear modes. The vanishing of A^^i^' at g > 1 
was explained analytically in our recent work [s^. In fact, the same explanation applies here, as the problem of the 
presence or absence of the finite threshold reduces to the existence or nonexistence of the solitons with a vanishingly 
small norm, N ^ 0. Obviously, in this limit the radial width of the mode trapped by the HO potential remains finite 
in the (x, y) plane, while the temporal width (in the r direction) becomes indefinitely large. Thus, the 3D problem 
for actually reduces to its 2D counterpart, for which an explanation of the absence of the threshold at g > 1 

was recently given in work [33| . 

We have also studied vortex solitons in the spatiotemporally anisotropic system, with Z? < 1 in Eqs. ([T]) (usually, 
the anomalous GVD is weaker at the SH [ij, therefore it makes sense to consider the case of D < 1). In particular. Fig. 
[3]Jc) shows, for this case, the vortex' slope at the center, defined in Eq. (fT5|) . as a function of A at g = 1 and 0.6 for 
the limit case oi D = (zero GVD at the SH component). Note that the vortex soliton disappears continuously both 
for q = 0.6 and 1 at Z? = [i.e., Nci{q — 0.6, D — 0) — 0], but it disappears by a jump at q — 0, i.e., Nci{q — 0, D ~ 0) 
is finite. 

The VA approximation poorly agrees with the numerical results for D ~ 0, which is explained by the fact that 
ansatz ([7]), postulating equal temporal widths of the FF and SH components of the vortex soliton, is not appropriate 
in this case. For D < 0, vortex-soliton solutions could not be found, in agreement with the known general principle 
that they do not exist in this case [1, |4lj . 



IV. STABILITY AND DYNAMICS OF THE VORTEX SOLITONS 



Systematic simulations demonstrate that the vortex solitons found above are stable in a large part of their existence 
domain. However, they are subject to an oscillatory instability at A^ > Nc2, when the total energy (norm) of the 
soliton exceeds the respective critical value, A"c2- An example of such an instability is presented in Fig. |31[a), which 




FIG. 3: (a) The slope of the vortical profile at the center, defined as per Eq. (|13p . produced by the numerical results (rhombuses), 
and by the VA (dashed curves) at q = 1 and 0.6 for D = 1. (b) The critical value, Nci, as predicted by the VA (the solid curve) 
and found from the numerical results (rhombuses), for D — 1. The vortex solitons do not exist at A'^ < Nc. (c) The same as in 
(a), but for D = (the case of zero GVD of the SH component). 




t t 



FIG. 4: (a) The perturbed evolution of the Fourier amplitudes iti and m_i, defined as per Eq. (|14[). for D = 1. (b) The same 
for D = 0. 



shows the evolution of two azimuthal Fourier amplitudes, which are defined as 



u±i{z) 



u{x, y, z, T)e^^^ dxdydr 



(14) 



dX q — l,D — \ and N — 185 (obviously, for the stationary vortex u+i = const, and m_i = 0). In the course of the 
unstable evolution, amplitude U-i{z) increases from zero, exhibiting oscillatory dynamics. Further, Fig. Sfb) shows 
the perturbed evolution of the same amplitudes iov q ^ 1, N ~ 185, and D = Q (zero GVD at the SH component). 
The difference from the picture shown inSJ^a) is that, at = 0, the oscillations are apparently chaotic. 

Figure [5] demonstrates that, in terms of the shape of the same vortex soliton whose evolution is illustrated by Fig. 
m^a), the oscillatory perturbations induce an azimuthal modulational instability of the axisymmetric shape. The axial 
symmetry is partly destroyed and then restored periodically, as one can see from the comparison of Figs. [5] (a) and 
(b). In the former panel, the modulation depth attains its maximum exactly at point z — 110, where U-i(z) has a 
maximum in Fig. |4] (a), and the axial symmetry is restored at z = 150, where it_i(z) vanishes in Fig. IHa). The 
period of the periodic breakup and retrieval of the axial symmetry is Z sa 80, the same as observed in the oscillations 




FIG. 5: Plots of \u{x,y) \ in the cross section drawn through r = T/2 at z = 110 (a) and z = 150 (b), for the oscillatory 
unstable vortex soliton with N = 185. 




FIG. 7; Collisions between moving vortex solitons with identical vorticities, under the action of kicks uj = ±0.5 (a); identical 
vorticities and u) = ±1 (b); opposite vorticities, ±1, and kicks u} — ±0.5 (c) and uj = ±1 (d). Shown are profiles in the cross 
section drawn in the longitudinal (r) direction through x — L/2 = 1.66 and y = L/2. 



of U-i{z) in Fig. Ilja). This dynamical regime resembles the ones featuring periodic splittings and recombinations of 
2D vortices, which are trapped in the isotropic HO potential in the systems vifith x*^^-' fssj and self-attractive x*^*^^ (29l | 
nonlinearities. 

The critical value for the onset of the oscillatory instability, found from the direct simulations, is shown in Fig. [nja) 
for several values of the mismatch, q = 0.2, 0.6, 1.0, 1.4, and 1.8, in the systems with the spatiotemporally isotropic 
diffraction-dispersion {D — 1), and with zero GVD of the SH component {D — 0). Stable vortex solitons exists in the 
broad region of Nci < N < Nc2 [recall Nci is the existence threshold for the solitons, see Fig. El^b)]. It is worthy to 
note that the stability region is somewhat larger for D = Q than for D = 1. 

The obvious Galilean invariance of Eqs. ([ij with D — 1 (the case of the spatiotemporally isotropic diffraction- 
dispersion operators) in the longitudinal (t) direction suggests a possibility to study collisions between moving solitons, 
if they are set in motion by means of the application of properly matched kicks to the FF and SH components, 
{u,v) (ue"*'^'^, ■ye"^*"'^), where w is an arbitrary boost parameter. Figure [Tj a) shows a typical example of the 
simulated head-on collision of two vortex solitons with equal vorticities 1, kicked by w = ±0.5. The collision is 
strongly inelastic, transforming the two colliding vortices into three, one quiescent and two symmetrically moving 
ones. We stress that each soliton produced by this inelastic interaction is a vortex with the same vorticity, 1, 
which was embedded into the initial modes, the total angular momentum ([3]) and the linear momenetum (j4]) being 
conserved. Inelastic collisions between vortices do not generate additional zero- vorticity solitons. The increase of the 
kick to Lo — ±1 results in a quasi-elastic passage of the colliding solitons through each other, see Fig. Eljb). Figures 
[TJc) and (d) demonstrate the collision of two solitons with opposite vorticities, ±1 and —1 kicked by (c) lo ~ ±0.5 
and (d) a; = ±l.@The collision is inelastic for w = ±0.5 and quasi-elastic for uj = ±1. 



V. CONCLUSION 



The objective of this work is to propose the model of optical media with the x^'^^ nonlinearity, in which 3D vortical 
solitons, which are always unstable against splitting in the free space, may be stabilized by means of the effective 2D 
isotropic HO (harmonic-oscillator) potential, induced by an appropriate spatial modulation of the refractive index. A 
broad stability region for the 3D vortices has been found, with their shape accurately predicted by the VA (variational 
approximation). The vortices exist with the norm (total energy) exceeding a finite threshold, which vanishes at values 



of the mismatch q > 1. The latter fact was explained by means of the comparison with similar results recently 
reported in the 2D version of the system in our work [33l |. The stability region is limited by another critical value 
of the norm, which is much higher than the existence threshold. Above the stability border, the vortex oscillates, 
periodically breaking and restoring its axisymmetric shape. Collisions between vortices, which may move freely in the 
longitudinal (temporal) direction, were also studied. If the collision velocity is small enough, the collision transforms 
the two colliding vortices into three, with the same vorticity. The results reported in this work suggest that 3D optical 
vortex solitons may be created in the stable form. 

The analysis can be naturally extended in other directions. In particular, it may be interesting to study the system 
with the Type-II x^^^ nonlinearity, i.e., the three- wave system with two FF components determined by their mutually 
orthogonal polarizations 0, Further, it is relevant to consider a possibility to stabilize the 3D vortex solitons using 
a periodic (lattice) potential, instead of the HO one, cf. the analysis performed previously for the x^^^ IHl and x^^^ 
[43l445l | nonlinearities. Also relevant may be the analysis of the stability of trapped vortices in a system including 
competing x'^-* and x^'^'' nonlinearities, cf. the consideration of the free-space model in work 
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